Statistics of Earthquakes in Simple Models of Heterogeneous Faults 
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Simple models for ruptures along a heterogeneous earthquake fault zone are studied, focussing on 
the interplay between the roles of disorder and dynamical effects. A class of models are found to 
operate naturally at a critical point whose properties yield power law scaling of earthquake statistics. 
Various dynamical effects can change the behavior to a distribution of small events combined with 
characteristic system size events. The studies employ various analytic methods as well as simulations. 



PACS numbers: 91.30.P, 62.20.F, 62.20.D 

The Gutenberg-Richter law [jD for the statistics of 
earthquakes — frequency inversely proportional to a 
power of the seismic moment — is well established over 
about 10 orders of magnitude. It is clearly a property of 
regional fault systems. The statistics of earthquakes on 
individual faults is much more controversial, indeed given 
the degree of geometrical complexity usually observed it 
is not even clear whether single faults are well defined. 
Nevertheless, statistics in various narrow fault zones in 
which slip is primarily along one direction — which we 
will henceforth refer to as "faults" — have been studied, 
and the behavior is found to vary substantially. In par- 
ticular, Wesnousky || has found that faults with large 
total displacement which are relatively regular typically 
have a power law distribution only for small events - 
if at all — and events with a much larger characteristic 
size in which the whole fault slips, with few events in be- 
tween. In contrast, less mature faults with more irregular 
geometries can have power law statistics over the whole 
range of observed magnitudes j2| . 

In this paper, we will show that simple models which 
include fault plane heterogeneities can exhibit both of 
these types of behavior and analyze the origin of the 
power law statistics and departures from it in these sys- 
tems. In particular, we will argue that power law statis- 
tics can be understood quantitatively in terms of proxim- 
ity to a specific non-equilibrium dynamical critical point. 
Like most critical points, the resulting exponents, al- 
though "universal" , will depend on certain properties in 
the system: the dimensionality, the range of interactions, 
randomness, and perhaps other aspects. 

Most previous work on simple models has involved 
variants of the Burridge-Knopoff (or "sliderblock" ) model 
in which the randomness is generated dynamically and 
inertia and friction laws play an essential role [^) . These 
systems appear to exhibit power law statistics over some 
range with a cutoff beyond some magnitude and with 
most of the slip occurring in larger, system size events. 
But the understanding of the origin of the power-law be- 
havior is very limited. Our approach here will be to start 



with analytic understanding of a class of models and then 
add in various additional physical features by analytic 
scaling arguments in the framework of the renormaliza- 
tion group (RG), aided by numerical studies. 

To investigate possible critical points, we first study 
infinite systems driven by a constant drive force F. The 
dynamical variables u(r, t) represent the discontinuity 
across the fault plane in the component of the displace- 
ment in the direction of slip. We consider general equa- 
tions of motion of the form: 

779u(r, t)/dt — F + a(r, t) - f R [u{r, t), r, {«(r, t' < t)}] 

(1) 

where 

t 

a(r,t)= J dt' J d d r'J(r-r',t-t') [u(r', t') - u{v,t)] 

- — OC 

(2) 

is the stress and fn is a quenched random "pinning" force 
crudely representing inhomogeneities in the friction, as- 
perities, stepovers etc., which in general can depend on 
the local past history (e.g. as in velocity dependent fric- 
tion). The dynamics will be determined by this local 
history dependence, the stress transfer function J(r,i), 
and the coefficient i] [II . 

Substantial simplifications occur if fn is history inde- 
pendent and J(r,i) > for all (r, t); we will call these 
monotonia models. Related monotonic models have been 
studied extensively in various other contexts m . Their 
crucial simplifying feature is that the steady state veloc- 
ity v = (du/dt) is a history independent function of F [Q. 
For F less than a critical force F c , v = 0, while just above 
F c , v ~ (F — F c )p . Universal scaling behavior exists on 
large length scales near F c . Quasi-static properties such 
as exponents and scaling functions depend on only a few 
quantities: the spatial dimension d; the range of the inter- 
actions if they are long-range, i.e. with the static stress 
transfer J s (r) = J dtJ(r,t) ~ l/r d+r , with T < 2; and 
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the range of correlations in fn, which we will generally 
assume are short-range in u and r. Long time dynamic 
properties such as f3 depend in addition on the small u> 
dependence of J(q, u). If F is adiabatically increased 
towards F c , the system moves from one metastable con- 
figuration to another by a sequence of "quakes" of various 
sizes. The "quakes" can be characterized by their radius 
R, the d-dimensional area A which slips (by more than 
some small cutoff), their moment M = J d d rAu(r), a 
typical displacement Au ~ M/A, and a duration r. 

From RG expansions B around a dynamic mean field 
version of Eq. ([!]) and scaling arguments it is found that 
for large quakes Au ~ R>, A ~ R di with df < d a fractal 
dimension, M ~ R d f+^ and r ~ R z . The distribution of 
moments is 

P(M)dM ~ dM / M 1+B p 00 (M / M) (3) 

with poo a universal scaling function which decays ex- 
ponentially for large argument. The cutoff M for large 
moments is characterized by a correlation length — the 
largest likely radius — £ ~ l/(F c -F) v with M ~ i df+C - 
In mean-field theory, B — 1/2, the quakes are fractal 
and displacements are of order the range of correlations 
in fn(u), i.e. £ = 0. The mean-field exponents are valid 
for d > 4(f) = 2f where f = min(T, 2) [§. For a planar 
fault in an elastic half space, d = 2 and r = 1 ||; the 
physical system is thus at the upper critical dimension 
d = d c §. 

As usual, at the upper critical dimension, there are log- 
arithmic corrections to mean-field results. We find barely 
fractal quakes with A ~ R 2 /\nR so that the fraction of 
the area slipped decreases only as 1/ In r away from the 
"hypocenter" . The typical slip is Au ~ (lni?) 1 / 3 so that 
M ~ i? 2 /(hii?) 2/3 - The scaling form of P(M) is the 
same as Eq. (§) with the mean-field p OQl although for 
M « M, P(M) ~ (InM) 1 / 3 /!/ 3 / 2 so that B will be 
virtually indistinguishable from 1/2. 

We now consider more realistic drive and finite-fault- 
size effects. Driving the fault by very slow motion far 
away from the fault is roughly equivalent to driving it 
with a weak spring, i.e. replacing F in Eq. (|l|) by 
F(r,t) = K[vt - u(r,t)}. With v -> the system 
must then operate with the spring stretched to make 
F(r, t) < F c at least on average; it will actually oper- 
ate just below F c . Under a small increase, AF, with 
constant force drive, (Au) « nAF J MP(M)dM with 
n the number of quakes per unit area per increase in 
AF; n{F) is non-singular at F c ||. The known scaling 
laws yield (Au) - AF£( 2t+ M 1 - B '> ~ AF£ for our case. 
For consistency, we must have in steady state with the 
spring drive, KvAt — AF — KAu so that the system 
will operate with a correlation length £ ~ l/K 1 ^, i.e. 
1 /K for our case. For a fault section with linear dimen- 
sions of order L, drive either from uniformly moving fault 
boundaries or from a distance ~ L perpendicularly away 



from the fault plane will be like K ~ 1/L so that the 
power-law quake distribution will extend out to roughly 
the system size £ ~ L. For smaller quakes, i.e. R << L, 
the behavior will be the same as in the infinite system 
with constant F drive, but the cutoff of the distribution 
of moments will be like Eq. (S) with a different cutoff 
function p that depends on the shape of the fault, how it 
is driven, and the boundary conditions. 

We have tested these conclusions numerically by study- 
ing a discrete space, time, and displacement version of a 
monotonic Eq. (|l|) with quasistatic stress transfer appro- 
priate to an elastic half space ||. The slip, u, is purely 
in the horizontal direction along the fault and fn[u(r)] 
is a series of equal height spikes with spacings which are 
a random function of r. When a(r,t) > fn[u(r,t], u(r) 
jumps to the next spike. The boundary conditions on the 
bottom and sides are uniform slip — (u = vt) with in- 
finitesimal v — and stress free on the top. The statistics 
of the moments of the quakes are shown by the triangles 
in Fig. |. 




moment (a.u.) 

FIG. 1. Histograms of moments for a simulation of a rect- 
angular fault with 32x128 cells for the discrete monotonic 
quasistatic model. Triangles: without dynamical weakening 
(e = 0). Diamonds: with dynamic weakening with e = 0.95. 
(e is defined in eq. (^).) The straight line indicates the pre- 
dicted slope B = 1/2. 

Although the uncertainties are appreciable, relatively 
good agreement is found with the prediction B = 1/2. A 
typical large quake is illustrated in Fig. ^| (a) ; it appears 
almost fractal as predicted and will tend to stay away 
from the bottom and sides. The ratios of the moments 
of quakes to their areas have been studied and found to 
grow only very slowly with the predicted. This is 

in striking contrast to earthquakes in conventional crack 
models which are compact and have Au ~ R (i.e. £ = 1), 
so that M I A ~ VA. 

Because the system is at its critical dimension, the cut- 
off function p of the moment distribution appropriate to 
the boundary conditions, as well as various aspects of the 
shapes and dynamics of quakes can be computed. For 
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quasistatic stress transfer, J(r, f) ~ 5(t)/r 3 , in the infi- 
nite system, the quake durations are found to be r ~ R z 
with z < 1 for d < d Cl corresponding to an unphysical su- 
personic propagation of disturbances || . In the marginal 
dimension d — d c , z = 1 with logarithmic corrections. A 
more physical dynamics with sound-travel-time delay has 
slower growth of the quakes with z = 1 in all dimensions. 
In either case the growth will be very irregular — in- 
cluding regions starting and stopping — in contrast to 
crack models and what is often assumed in seismological 
analyses of earthquakes. 




FIG. 2. Distribution of horizontal slip, u, along a fault with 
32x128 cells for a single large quake event. The lighter a cell 
the bigger its slip during the quake. Top (a) almost fractal 
quake with a total moment of 1750 (and 1691 cells failing) for 
the monotonic model without any dynamical effects (e = 0). 
Bottom (b) "cracklike" quake with a total moment of 16922 
(and 2095 cells failing) for the model with dynamic weakening 
(e = 0.95). In both cases the system is driven by horizontally 
creeping fault boundaries (sides and bottom) while the top 
boundary is free. 

A crucial feature of monotonic models is that the slip 
profile Att(r) of a quake is independent of the dynamics 
0. But the most interesting dynamical issues concern 
the effects left out of the monotonic models that can make 
this feature breakdown. We first consider including some 
weakening effects of sections which have already slipped 
in a given quake. This is best studied in the discrete 
model. To crudely model a difference between static and 
dynamic friction we choose, 

f R = /a[u(r),r]{l - eQ[u(r, t) - u(r,t - T)}} (4) 

with T a cutoff time much longer than the duration of 
the largest quakes, but much smaller than the interval 
between the quakes. This can crudely represent a differ- 
ence between static and dynamic friction. For e = and 
non-negative J the model is monotonic, while for e > 
it is non-monotonic. The effects of small weakening can 
be analyzed perturbatively. 

With e = 0, consider a quake of diameter Ri (<< L 
or £), with moment Mi and area A±: i.e. A\ sites have 
slipped. If a small e is turned on at the end of the quake, 
all slipped sites that are within e of slipping will now 
slip again — this will be N^" ~ tA\ sites. The sim- 
plest justifiable guess is that each of these will cause 



an approximately independent secondary quake. The 
total moment of these secondary quakes will be domi- 
nated by the largest one, so the extra moment will be 
M% x ~ (eA-t) 1 /®. If Mf* << Mi this process can con- 
tinue but will not increase the total moment substan- 
tially. But if M% x ~ Mi, the process can continue with 
a larger area A^ and hence a larger M ex , leading to run- 
away. The scaling laws yield B = (df —T + ()/(df + () 
and A ~ M d f^ d f + ^ with C < f , so that for any e, 
for large enough Mi, Ml > M D ~ e -( rf /+0/(f-C) ; m%* 
will be comparable to Mi and the quake will become 
much larger. In the case of interest Mo ~ e~ 2 . In the 
force driven infinite system for F < F c , quakes of size £ 
will runaway and become infinite if £ > e -1 '^^^'. Since 
f ~ (F - F c )~ v and 1/v = f - £ this will occur for 
F c — F < C w e with some constant C w . This result is very 
intuitive and justifies a posteriori the assumptions lead- 
ing to it: Since on slipping, the random pinning forces, 
fn in a region are reduced by order e, the effective criti- 
cal force F c for continuous slip will have been reduced by 
order e; thus if F > F c (e) = F c — C w e, the mean velocity 
v will be nonzero. 

A similar but more subtle effect can be caused by stress 
pulses that result from non-positive J(r,t); these arise 
naturally when one includes elastodynamic effects. We 
consider 

J(r,i) ~ 5(t-r/c)/r d+r + aS\t-r/c)/cr d+ ~< (5) 

with c the sound speed. The scalar approximation to 
elasticity in a half space corresponds to d = 2, T = 1, 
7 = 0, and a = 1. If a region slips forward, the stress at 
another point first has a short pulse at the sound arrival 
time from the second term in Eq. (^|), then settles down 
to its smaller static value, i.e. it is non-monotonic. The 
magnitude of these stress pulses and their duration is set 
by various aspects of the models, for example larger rj in 
Eq. |l]) implies weaker stress pulses as the local motion 
will be slower. By considering which of the sites in a 
long quake with a = can be caused to slip further 
by such stress pulses — here the dynamics matters - 
we find that runaway will occur for M > Mjj ~ a~ 4 
for the physical case. We have checked this in d = 1 
with r = 1 and 7 = finding the predicted reduced 
critical force F c (a) ~ F c — C v o? as shown in Fig. ||. 
These 1-d simulations also reveal a hysteretic v(F) curve 
in finite systems. This should also occur with the velocity 
weakening model. 

We can now understand what should happen with ei- 
ther weakening or stress pulses in finite systems driven 
with a weak spring or with slowly moving boundaries. 
As the system is loaded, quakes of increasing size can 
occur. If the system is small enough that it cannot 
sustain quakes with M > M£>(e,a) then the behavior 
will not be much different from the monotonic case with 
e = a = 0. This will occur if L < R D (e,a) - M D /2 ~ 
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max(C a /a 2 ,C e /e) with appropriate coefficients C a , C e , 
which will depend on the amount of randomness in the 
fault. But if L > Ftp, quakes of size of order Ro will 
runaway and most of the system will slip, stopping only 
when the load has decreased enough to make the loading 
forces less than the lower end of the hysteresis loop in 
v{F) (as in Fig. |). 
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be understood by scaling arguments. It is hoped that 
other physical phenomena such as geometrical disorder, 
side branching, and multiple cracks might start to be 
addressed in this framework. We note one extra effect 
which can be readily analyzed: long range correlations in 
the randomness (perhaps caused by prior history of the 
fault) . Varying the power-law of the decay of correlations 
of fn increases B continually from 1/2 to 2/3 and £ con- 
comitantly from to 1 with the quakes becoming more 
compact and crack-like as the randomness correlations 
become longer range. 
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FIG. 3. Mean velocity versus force for one dimensional sys- 
tem with a non-monotonic kernel 
J{x,t) = S(t - x)/x 2 + aS'(t - as) /as for a = 0.8,0.5,0. A 
spring or boundary loaded system will traverse the hysteresis 
loops in the direction indicated. Inset: the threshold force, 
on increasing the load; A a = [l~Fj (a)/ Fj (a = 0)] 1/2 
is plotted vs. a. 

Because of the tendency of regions that have already 
slipped to slip further, and the consequent buildup of 
larger stresses near the boundaries of the slipped regions, 
large events in systems with dynamic weakening will be 
much more cracklike than in monotonic models, probably 
with Au ~ L. Statistics of quakes with weakening, e, 
reasonably large, but no stress pulses (a — 0) are shown 
in Fig. [j] and in ||; note the absence of quakes with 
intermediate moments. A typical large event in this case 
is shown in Fig. || (b); it appears to be crack-like. 

In this paper we have shown that simple models of 
heterogeneous faults — with the dimensionality and long- 
range elastic interactions properly included — can give 
rise to either power-law statistics of earthquake moments 
or a distribution of small events combined with charac- 
teristic system size events. Which behavior — or in- 
termediate behavior — obtains is found to depend on a 
number of physical properties such as frictional weaken- 
ing and dynamic stress transfer, analogs of which should 
definitely exist in real systems. In the power-law-regime 
the conventionally defined Gutenberg-Richter exponent 
b = 3.B/2 is is found to be b — 3/4. This is intriguingly 
close to values observed by Wesnousky ||, but it is not 
clear if any significance should be attached to this. 

More significant is the framework that we have built, 
which enables certain results (and many more not pre- 
sented here) to be obtained analytically and others to 
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